
net get eacsap.pkg, from(http://www.stata.com/data/jwooldridge/)
use airfare.dta, clear
egen double concenbar = mean(concen), by(id)
sum concenbar
gen double concenbar_dm = concenbar - r(mean)
sum ldist
gen double ldist_dm = ldist - r(mean)

xtreg lfare concen y98 y99 y00, fe vce(cluster id)

* Usual Mundlak regression, with time constant variables added:
reg lfare concen concenbar ldist c.ldist_dm#c.ldist_dm ///
    y98 y99 y00, vce(cluster id)     
		

* Heterogeneous time effects using Mundlak:
reg lfare concen concenbar c.ldist c.ldist_dm#c.ldist_dm  /// 
    y98 y99 y00 c.y98#c.concenbar_dm c.y99#c.concenbar_dm ///
		c.y00#c.concenbar_dm, vce(cluster id)
	

* FE gives same estimates:
xtreg lfare concen y98 y99 y00 c.y98#c.concenbar_dm ///
  c.y99#c.concenbar_dm c.y00#c.concenbar_dm, fe vce(cluster id)
	
	
* Interact distance with time dummies:
reg lfare concen concenbar ldist c.ldist_dm#c.ldist_dm     /// 
    y98 y99 y00 c.y98#c.concenbar_dm c.y99#c.concenbar_dm  ///
    c.y00#c.concenbar_dm c.y98#c.ldist_dm c.y99#c.ldist_dm ///
    c.y00#c.ldist_dm, vce(cluster id)
		
test c.y98#c.concenbar_dm c.y99#c.concenbar_dm ///
     c.y00#c.concenbar_dm c.y98#c.ldist_dm     ///
		 c.y99#c.ldist_dm c.y00#c.ldist_dm

		
* Now heterogeneous slope on concen, too.
reg lfare concen c.concen#c.concenbar_dm c.concen#c.ldist_dm ///
    concenbar ldist c.ldist_dm#c.ldist_dm y98 y99 y00 ///
    c.concenbar#c.concenbar_dm c.concenbar#c.ldist_dm ///
    c.y98#c.concenbar_dm c.y99#c.concenbar_dm ///
    c.y00#c.concenbar_dm c.y98#c.ldist_dm     ///
    c.y99#c.ldist_dm c.y00#c.ldist_dm, vce(cluster id)
	

egen lfarebar = mean(lfare), by(id)
sum lfarebar
gen lfarebar_dm = lfarebar - r(mean)
* First apply the previous model where only time averages are used.
reg lpassen concen concenbar ldist ldistsq y98 y99   ///
    y00 c.y98#c.concenbar_dm c.y99#c.concenbar_dm    ///
    c.y00#c.concenbar_dm lfare lfarebar c.y98#c.lfarebar_dm /// 
    c.y99#c.lfarebar_dm c.y00#c.lfarebar_dm, vce(cluster id)
	
* Now use the unit-specific linear trend heterogeneity:	
gen double lpassen_h = .
gen double lpassen_dt = .
gen double lfare_h = .
gen double lfare_dt = .
gen double lfare_a0 = .
gen double lfare_a1 = .
gen double concen_h = .
gen double concen_dt = .
gen double concen_a0 = .
gen double concen_a1 = .
gen double y98_h = .
gen double y99_h = .
gen double y00_h = .
gen double y98_dt = .
gen double y99_dt = .
gen double y00_dt = .
gen double y98_a1_t = .
gen double y99_a1_t = .
gen double y00_a1_t = .

gen t = year - 1997

local i = 1

while `i' <= 1149 {

	qui reg lpassen t if id == `i'
		predict lpassen_t, xb
		qui replace lpassen_h =lpassen_t if id == `i'
		qui replace lpassen_dt =lpassen - lpassen_t if id == `i'
		
	qui reg lfare t if id== `i'
		qui replace lfare_a0 =  _b[_cons] if id ==  `i'
		qui replace lfare_a1 = _b[t] if id == `i'
		predict lfare_t, xb
		qui replace lfare_h = lfare_t if id == `i'
		qui replace lfare_dt = lfare - lfare_t if id == `i'
		
	qui reg concen t if id == `i'
		qui replace concen_a0 = _b[_cons] if id == `i'
		qui replace concen_a1 = _b[t] if id == `i'
		predict concen_t, xb
		qui replace concen_h = concen_t if id == `i'
	
	qui replace concen_dt = concen - concen_t if id == `i'
		qui reg y98 t if id == `i'
		predict y98_t, xb
			qui replace y98_h = y98_t if id ==`i'
		qui replace y98_dt = y98 - y98_t if id == `i'
			qui replace y98_a1_t = _b[t]*t if id == `i'
			
		 qui reg y99 t if id == `i'
		predict y99_t,xb
			qui replace y99_h = y99_t if id == `i'
		qui replace y99_dt = y99 - y99_t if id == `i'
			qui replace y99_a1_t = _b[t]*t if id == `i'
		
		qui reg y00 t if id == `i'
		predict y00_t, xb
			qui replace y00_h = y00_t if id == `i'
		qui replace y00_dt = y00 - y00_t if id == `i'
			qui replace y00_a1_t = _b[t]*t if id == `i'
		
		drop lfare_t concen_t lpassen_t y98_t y99_t y00_t
		
		local i = `i' + 1
		}
	

sum concen_h
gen double concen_h_dm = concen_h - r(mean)

sum lfare_h
gen double lfare_h_dm = lfare_h - r(mean)

sum concen_a0
gen double concen_a0_dm = concen_a0 - r(mean)

sum concen_a1
gen double concen_a1_dm = concen_a1 - r(mean)

sum lfare_a0
gen double lfare_a0_dm = lfare_a0 - r(mean)

sum lfare_a1
gen double lfare_a1_dm = lfare_a1 - r(mean)
		
		
reg lpassen_dt concen_dt lfare_dt y98_dt y99_dt y00_dt ///
    c.y98_dt#c.concen_a0_dm c.y99_dt#c.concen_a0_dm    ///
    c.y00_dt#c.concen_a0_dm c.y98_dt#c.lfare_a0_dm     ///
    c.y99_dt#c.lfare_a0_dm c.y00_dt#c.lfare_a0_dm      ///
    c.y98_dt#c.concen_a1_dm c.y99_dt#c.concen_a1_dm    ///
    c.y00_dt#c.concen_a1_dm c.y98_dt#c.lfare_a1_dm     ///
    c.y99_dt#c.lfare_a1_dm c.y00_dt#c.lfare_a1_dm,     ///
    nocons vce(cluster id)
	
reg lpassen concen lfare c.ldist c.ldist_dm#c.ldist_dm  ///
    y98 y99 y00 concenbar lfarebar concen_h lfare_h     ///
    c.y98#c.concen_a0_dm c.y99#c.concen_a0_dm           ///
    c.y00#c.concen_a0_dm c.y98#c.lfare_a0_dm            ///
    c.y99#c.lfare_a0_dm c.y00#c.lfare_a0_dm             ///
    c.y98#c.concen_a1_dm c.y99#c.concen_a1_dm           ///
    c.y00#c.concen_a1_dm c.y98#c.lfare_a1_dm            ///
    c.y99#c.lfare_a1_dm c.y00#c.lfare_a1_dm,            ///
    vce(cluster id)
	
	
